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Three programs in Mathematica are presented, which produce expressions for the lowest order 
and the higher order corrections of the Phase Integral Approximation. First program is pertinent 
to one ordinary differential equation of the Schrodinger type. The remaining two refer to a set of 
two such equations. 



I. INTRODUCTION 

In this paper we present three programs in Mathematica related to the Phase Integral Approximation (PIA) [H. 
They produce expressions for the lowest order aproximation and the higher order corrections. The first program 
generates the higher order corrections Y2n{x) pertinent to one ordinary differential equation of the Schrodinger type. 
The second program gives the vectors hm{x) and the third one the corrections Ym{x) and the vectors Srn{x). These 
quantities are pertinent to a set of two ordinary differential equation of the Schrodinger type. The programs will 
be referred to as Y2n, bvm and Ymsvm, respectively. They should be saved in file. mat (file equal to Y2n, bvm 
or Ymsvm) and run by using the standard Mathematica input command << file. mat. The user supplied input data 
should be saved in file . dat. If file . dat is an empty file, the default input data contained in each program will be 
used in computation. 

Each program opens a dialog in Mathematica session in which the user is asked to specify the form of output from 
Mathematica (DutputForm, TeXForm or FortranForm), give the number of terms to be determined and answer a few 
questions specific to each program. The results produced by Mathematica will be saved as file. res, file.resTeX 
or file.resFor. All equation and section numbers given in what follows refer to [1]. 

One should answer y to the dialog question Expand [ , Trig -> True ] , y/n?, if the user supplied data file 
contains some trig functions, and n otherwise. 

The factor g{x) after the dialog question Factor g[x] in eigenvector = ? in the program Ymsvm, can be chosen 
either to be equal to one, or to the denominator in the just printed expression for s02/s01 (= smix) / sm{x)) , see 
Eq. (124). In the last case, the implied factors should also be included, e.g., smx comming from tanx etc. 

All three programs deal with multiple sums, see Eqs. (40) and (54). These sums are programmed in the simplest 
possible way, e.g., the sum 

a+/3+7+(5+cr— m 

present in Eq. (54), where < a, /?, 7, (5, cr < m — 1, is programmed as 
sum2 = ; 

Do[ sum2 += If [ a + b + g + d + s == m, (*then*) 
Y[a] Y[b] Y[g] Y[d] sv[s], (*else*) ], 

{a, 0, m - 1}, {b, 0, m - 1}, {g, 0, m - 1}, {d, 0, m - 1}, {s, 1, m - 1} ]; 

etc. This makes the programs as close as possible to mathematical formulas thereby eliminating programming errors. 
For the same reasons, the integral which defines the coordinate (e,Sm), see Eq. (105), was not simplified by using 
Eqs. (106)-(108) which would make the computation faster. However, this would make programming a bit more 
complicated and error-prone. In our computations, this type od optimalizations was not necessary. Our aim was to 
produce correct results in a reasonable time (seconds or minutes rather than hours). A strong test for the correctness 
of programming was the vanishing of the odd order corrections, Y2n-i{x) = (which required cancellation of many 
terms), see Sec. VIIIA. Another check was the fact that in hermitian cases, the same results were produced in the 
simplified hermitian and non-hermitian theory, see Sees. VIIIA- VIIIC and VIIIE. 
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II. PROGRAM TO DETERMINE THE CORECTIONS Yzn FROM THE RECURRENCE RELATION (40) 
File Y2n.inat: 

Calculation of the phase integral corrections Y2n from recurrence relations, 
Eq. (40) in [1], for a scalar case, i.e. for one ODE of the Schroedinger type: 

u' ' (x) + R(x) u(x) = 0. 

[1] A. A. Skorupski, "Phase Integral Approximation for coupled ODEs of the 
Schroedinger type" , arXiv: 0710.5868, Sec. II. 

***************** Define type of output from Mathematica ********************* 
*) 

outpform = InputString ["Output, TeX or Fortrain form of results, o/t/f? "] ; 
sc = If[ outpform == "o", DpenWrite ["Y2n . res" , FormatType -> DutputForm] , 
If [ outpform == "t", OpenWrite ["Y2n.resTeX" , FormatType -> TeXForm] , 
OpenWrite ["Y2n.resFor" , FormatType -> FortranForm ] ] ]; 

(**) 

outpYm = InputString [ 

"Simple fractions or Common denominator or in Y2n, s/c? "] ; 

(**) 

WriteString [sc , "\n Formulas for corrections Y2n as functions of x or z 

(= zeta variable) . \n"] ; 

(* 

****************** Define maximum value of n in Y[2 n] ********************* 

*) 

nmax = Input ["nmax = ? "] ; 

WriteString [sc, "\n nmax = "] ; Write [ sc, nmax ]; 
(* 

***************** Define type of input to Mathematica ********************* 
*) 

inptform = InputString ["Input of R(x) and a(x) or General Y2n, i/g? "] ; 
(**) 

to = TimeUsed[] ; 
(**) 

If [ inptform == "i", 
(*then*) 
(* 

********************** Define default input data *************************** 
*) 

(*** Parabolic Model ***) 
af = 0; 

R = coef (x~2 - xl"2) ; 
(* 

**************** Read new input data from file Y2n.dat ********************* 
*) 

« Y2n.dat; 

WriteString [sc, "\n Y2n[x] for \n"] ; 
(* 

*) 

WriteString [sc, "\n R[x] = \n"] ; Write[sc, R] ; 

WriteString [sc, "\n Auxiliary function a[x] = \n"] ; Write[sc, af ] ; 



3 

(**) 

Qsq = R - af ; dQsq = D[ Qsq, x ] ; 
Qsqorl = Qsq; 

epO = ( (5/16) (dQsq/Qsq)~2 - (1/4) D[ dQsq, x ]/Qsq + af )/Qsq; 
aux = Simplify [ epO ]; 
auxl = Together [aux] ; 

WriteString[sc, "\n epsO[x] = \n"] ; Write[sc, auxl], 

(*else*) 

(********** Prepare quantities for general calculations 
(**) 

epO = epsO[x]; Qsq = Qsqr[x]; 

xorz = InputString ["x or zeta variable, x/z? "] ; 
If [ xorz == "x", (*then*) WriteString [sc , 

"\n Y2n[x] as functions of epsO[x], Qsqr [x] = Q"2[x] and derivatives \n"] ; 

Qsqorl = Qsq, (*else*) 

Qsqorl = 1; x = z; WriteString [sc , 

"\n Y2n[z] as functions of epsO[z] and derivatives \n"] ; ] 

]; 

(**) 

Qm2 = 1/Qsqorl; 
Y[0] = 1; 

(* 

********************** Start iterations for Y[2 n] ************************* 

*) 

For[n = 1, n <= nmeix, n++, 

suml = 0; sum2 = 0; sum3 = 0; 
m = 2 n; 

Do[ suml += If [ a + b == m, (*then*) Y[a] Y[b], (*else*) ], 

{a, 0, m - 2, 2>, -[b, 0, m - 2, 2} ] ; 
Do[ sum2 += If [ a + b + g + d == m, (*then*) 

Y[a] Y[b] Y[g] Y[d], (*else*) ], 

{a, 0, m - 2, 2>, ih, 0, m - 2, 2}, {g, 0, m - 2, 2}, {d, 0, m - 2, 2} ] ; 

Do[ sum3 += If [ a + b == m - 2, (*then*) 

epO Y[a] Y[b] + (3/4) Qm2 D[Y[a], x] D[Y[b], x] - 

(1/2) Y[a] qm2 (D[Y[b], {x, 2>] - (1/2) Qm2 D[ Qsqorl, x] D[Y[b], x] ), 
(*else*) ], {a, 0, m - 2, 2}, {b, 0, m - 2, 2} ] ; 

(**) 

Y[m] = (1/2) (suml - sum2 + sum3) ; 
(**) 

]; 

t = TimeUsed [] ; 

WriteString [sc, "\n CPU time used for computation (seconds) = "] ; 

Write [sc, t - tO] ; 

(* 

*********************** Simplify and write results *************************** 

*) 

For [n = 1 , n <= nmax , n++ , 
m = 2 n; 

WriteString [sc, "\n n = "] ; Write[sc, n] ; 
aux = Simplify [ Y[m] ] ; 

auxl = If[ outpYm == "c". Together [aux] , Apart [aux, x] ]; 
WriteString [sc, "\n Y2n = \n"] ; Write[ sc, auxl ] 

]; 

t = TimeUsed [] ; 

WriteString [sc , 

"\n CPU time used for computation & simplification (seconds) = "] ; 
Write [sc, t - tO] ; 
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End of file Y2n.mat. 

The file that follows is an example of the data file for the program Y2n. As it stands it is an ampty file containing 
only commemts. By uncommenting the definition of the functions a{x) and R{x): (* — > (**) and *) — > (**), one 
activates the input data. 

File Y2n.dat: 

(* Data pertinent to R[x] in the differential equation u' ' [x] + R[x] u[x] = 0. 

By default the auxiliary function af [x] = 0. For other choice 
include the data command: af = your_f unction [x] ; 

*) 

(************************ Budden'S Model: !|!!|!!|!!|!!|!!|!!|!!|!!|!!|!!|! ******* 

(* 

af = 0; 

R = coef x/(x - p) ; 
*) 

End of file Y2n.dat. 



III. PROGRAM TO DETERMINE FROM THE RECURRENCE RELATION (54) 

File bvm.mat: 

Calculation of bv[m] from the recurrence relation, Eq. (54) in [1] 

[1] A. A. Skorupski, "Phase Integral Approximation for coupled DDEs of the 
Schroedinger type", arXiv: 0710.5868, Sec. III. 

*) 

bvm = OpenWrite ["bvm.res" , FormatType -> OutputForm] ; 

(**) 

WriteString [bvm, "\n Formulas for vectors bvm as functions of x or z 

(= zeta variable) . \n"] ; 

to = TimeUsed[] ; 

Unprotect [Sqrt , Power] ; 

Sqrt[x_~2] := x; 

i"n_ := (-l)"(n/2) /; EvenQ [n] ; 

i-n_ := i (-l)-((n-l)/2) /; OddQ [n] ; 

re [x_] := Coef f icient [x, i, 0]; 

im[x_] := Coef f icient [x, i, 1]; 

cc[x_] := re[x] - i im[x] ; 

Potect [Sqrt, Power]; 

(* 

******************** Define maximum value of m in bv[m] ******************** 
*) 

mmax = Input ["\n mmax = ? "] ; 
(**) 

xorz = InputString["\n x or zeta variable, x/z? "] ; 
If[ xorz == "z" , (*then*) Q[x_] := 1; x = z ]; 
Qml = C![x]"(-1); Qm2 = C!ml"2; dQ = D [ Q [x] , x ] ; 
(**) 

Y[x, 0] = 1; 

bv[l] = i Qml D[sv[x, 0], x] ; 

YleqOQ = InputString ["\n Y[x, 1] = 0, y/n? "] ; 
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If[ YleqOQ == "y", (*then*) Y[x, 1] = ]; 

(* 

************************ Start iterations for bv[iii] ************************ 
*) 

For [m = 2 , m <= mmax + 1 , m++ , 
(**) 

siiml = 0; sum2 = 0; sum3 = 0; sum4 = 0; sum5 = 0; sum6 = 0; 
Do[ suml += If [ a + b + s == m, (*then*) 

Y[x, a] Y[x, b] ( sv[x, s] + 2 (Y[x, s] sv[x, 0] - bv[s]) ), (*else*) 

], {a, 0, m - 1}, -Cb, 0, m - 1}, {s, 1, m - 1} ]; 
Do[ sum2 += If [ a + b + g + d + s == m, (*tlien*) 

Y[x, a] Y[x, b] Y[x, g] Y[x, d] sv[x, s] , (*else*) ], 

{a, 0, m - 1}, {b, 0, m - 1}, {g, 0, m - 1}, {d, 0, m - 1}, {s, 1, m - 1} ]; 
Do[ sum3 += If[ a + b == m, (*then*) Y[x, a] Y[x, b] , (*else*) ], 

{a, 1, m - 1}, {b, 1, m - 1} ]; 
Do[ suin4 += If [ a + b + g + d == m, (*then*) 
Y[x, a] Y[x, b] Y[x, g] Y[x, d] , (*else*) ], 
{a, 0, m - 1}, {b, 0, m - 1>, {g, 0, m - 1}, {d, 0, m - 1} ] ; 
Do[ sums +=If[a + b + g + s==in-l, (*theii*) 

Y[x, a] Y[x, b] Y[x, g] Qml D[sv[x, s] , x] , (*else*) ], 
{a, 0, m - 1}, {b, 0, m - 1}, {g, 0, m - 1}, {s, 0, m - 1} ] ; 
Do[ sum6 += If [ a + b + s == m - 2, (*then*) 

Y[x, a] ( Y[x, b] ( qm2 ( D[sv[x, s] , {x, 2}] - Qml dQ D[sv[x, s] , x ] ) + 
epsO[x] sv[x, s] ) - Qm2 D[Y[x, b] , x] D[sv[x, s] , x] - 

(1/2) Qm2 ( D[Y[x, b] , {x, 2}] - Qml dQ D[Y[x, b] , x ] ) sv[x, s] ) + 
(3/4) Qm2 D[Y[x, a], x] D[Y[x, b] , x] sv[x, s] , (*else*) ], 
{a, 0, m - 2}, {b, 0, m - 2>, {s, 0, m - 2} ] ; 
(**) 

bv[m] = (1/2) (suml - sum2 + ( sum3 - sum4 ) sv[x, 0] + i 2 sum5 + sum6) ; 
(**) 

]; 

t = TimeUsed [] ; 

WriteString[bvm, "\n CPU time used (seconds) = "] ; Write[bvm, t - tO] ; 

(**) 

(* 

*********************** Simplify and write results *************************** 

*) 

For [m = 1 , m <= mmax , m++ , 

WriteString[bvm, "\n m = "] ; Write[bvm, m] ; 
aux = Simplify [ bv[m] ]; 
auxl = Together [aux] ; 

WriteStringCbvm, "\n bvm = \n"] ; Write [ bvm, auxl ]; 
(**) 

]; 

End of file bvm. mat. 



IV. PROGRAM TO DETERMINE Ym AND s„ FROM RECURRENCE RELATIONS FOR 2 ODES 
WITH EITHER HERMITIAN OR NON-HERMITIAN MATRIX, SEE [1], SEC. VI 

File Ymsvm.mat: 

Calculation of Y[m] and sv[m] from recurrence relations 
for 2 DDEs with either hermitiein or non-hermitiein matrix [1] . 



[1] A. A. Skorupski , "Phase Integral Approximation for coupled DDEs of the 
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Schroedinger type", arXiv: 0710.5868, Sec. VI. 

*) 

outpform = InputString ["Output, TeX or Fortrein form of results, o/t/f? "] ; 
vc = If [ outpform == "o", (*then*) OpenWrite ["Ymsvm.res" , 
Format Type -> OutputForm] , 

(*else*) If[ outpform == "t", (*then*) OpenWrite ["Ymsvm.resTeX" , 

FormatType -> TeXForm] , (*else*) OpenWrite ["Ymsvm.resFor" , 

Format Type -> FortrsinForm] ] ] ; 

to = TimeUsed [] ; 

Unprotect [Sqrt , Power] ; 

Sqrt[x_'2] := x; 

i"n_ := (-l)~(n/2) /; EvenQ [n] ; 

i-n_ := i (-l)-((n-l)/2) /; OddQ [n] ; 

re [x_] := Coef f icient [x, i, 0]; 

im[x_] := Coef f icient [x, i, 1]; 

cc [x_] := re [x] - i im[x] ; 

Potect [Sqrt , Power] ; 

(* 

*************** Define maximum value of m in Y[m] and sv[m] ****************** 
*) 

mmax = Input ["mmax = ? "] ; 

WriteString [vc , "\n mmax = "] ; Write[vc, mmax]; 

mmxpl = mmax + 1; 

(* 

*********************** Define default input data **************************** 

One must define the auxiliary function a(x, p) -> af eind the matrix elements 
Rjk(x, p) -> Rjk, where p represents parameter (s) . By default, the variable 
automatic = True. In that case, one must define a list of meaningful numerical 
replacements: parrepls = { x -> xO, p -> pO } which is necessary to calculate 
automatically the variables: Delta given by Eq. (121), sqrtDel (= Sqrt [Delta] ) 
and signQsq (= sign of Qsq given by Eq. (121)). If one puts automatic = False 
in the data file Ymsvm.dat, the definitions of Delta, sqrtDel and signQsq must 
be given in the file Ymsvm.dat. 
*) 

(*** ExEimple A ***) 
af = 0; 

Rll = x Cos[x]~2 + Sin[x]'"2; 
R12 = (x - 1) Cos[x] Sin[x]; 
R21 = (x - 1) Cos[x] Sin[x]; 
R22 = X Sin[x]~2 + Cos[x]"2; 
parrepls = { x -> 2 }; 
automatic = True; 
(* 

***************** Read new input data from file Ymsvm.dat ******************** 
*) 

<< Ymsvm.dat; 
(* 

*) ^ 



WriteString [vc , 


"\n 


Rll = 


\n"] ; 


Write [vc , 


Rll] 


WriteString [vc , 


"\n 


R12 = 


\n"] ; 


Write [vc. 


R12] 


WriteString [vc , 


"\n R21 = 


\n"] ; 


Write [vc. 


R21] 


WriteString [vc , 


"\n 


R22 = 


\n"] ; 


Write [vc , 


R22] 


WriteString [vc , 


"\n 


af = 


\n"] ; 


Write [vc. 


af] 



If [ automatic, (*then*) 
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WriteString[vc, "\n paramter replacement list = \n"] ; Write[vc, parrepls] , 
(*else*) 

WriteString [vc , "\n *** Non-automatic calculation *** \n"] 

]; 
(**) 

Gil = Rll - af ; G12 = R12; 021 = R21; G22 = R22 - af ; 

(**) 

exptrig = InputString ["Expand [ , Trig -> True ] , y/n? "] ; 

QsqmQ = InputString ["Qsqr with minus or plus Sqrt [ Delta ], m/p? "] ; 

WriteString [vc, "\n Qsqr with "] ; 

If [ QsqmQ == "m" , WriteString [vc , "minus "] , WriteString [vc , "plus "] ]; 

WriteString [vc, "Sqrt[ Delta] "] ; 

(**) 

(* 

************ Find eigenvalues and eigenvectors of the G matrix *************** 
*) 

If[ automatic, (*then*) Delta = (Gil - G22)"2 + 4 G12 G21 ]; 
If[ exptrig == "y". Delta = Expand[ Delta, Trig -> True ] ]; 
Delta = Factor [ Delta ] ; 

WriteString [vc, "\n Delta = \n"] ; Write[vc, Delta]; 
If [ automatic, (*then*) sqrtDel = Sqrt [ Delta ] ]; 
(**) 

Qsq = (1/2) (Gil + G22 + If [ QsqmQ == "m" , (*then*) -sqrtDel, (*else*) 
sqrtDel ] ) ; 

Qsq = Simplify [ Qsq ] ; 

If [ exptrig == "y", Qsq = Expand [ Qsq, Trig -> True ] ]; 
(**) 

If [ automatic, (*then*) 

signQsq = If [ ( Qsq /. parrepls ) < 0, (*then*) -1, (*else*) 1, 

(*aiid if neither True or False then*) 1 ] ] ; 
WriteString [vc , "\n signQsq = "] ; Write[vc, signQsq]; 
WriteString [vc, "\n Qsq = \n"] ; Write[vc, Qsq]; 
(**) 

(*** epsO = Simplify [ (Qsq" (1/4) D [Qsq" (-1/4) , {x, 2}] + af)/Qsq ]; ***) 
epsO = Together [ Simplify [ ( (5/16) ( D[qsq, x]/Qsq)~2 - 

(1/4) D[Qsq, {x, 2>]/Qsq + af )/Qsq ] ]; 
WriteString [vc, "\n epsO = \n"] ; Write[vc, epsO] ; 
(* « epsO.mat; *) 
(**) 

Q = If [ SignQsq < 0, (*then*) - i (- Qsq) "(1/2), (*else*) Qsq" (1/2) ]; 
WriteString [vc, "\n Q = \n"] ; Write[vc, Q] ; 

(**) 

Qml = Q"(-l); Qm2 = Qml"2; dQ = D[ Q, x ]; 
(**) 

s02os01 = ( Qsq - Gil )/G12; 

If[ exptrig == "y", s02os01 = Expand [ s02os01. Trig -> True ] ]; 
s02os01 = Simplify [ s02os01 ]; 

WriteString [vc, "\n s02/s01 = \n"] ; Write[vc, s02os01] ; 
Print [ "s02/s01 = ", s02os01 ]; 

(**) 

fact = Input ["Factor g[x] in eigenvector = ? "] ; 
WriteString [vc, "\n Factor g[x] = \n"] ; Write[vc, fact]; 
sOvl = fact; s0v2 = fact s02os01; 
(**) 

asqr = sOvl cc [sOvl] + s0v2 cc [s0v2] ; 

If [ exptrig == "y", asqr = Expand[ asqr. Trig -> True ] ] ; 
(**) 

normeigv = InputString ["Normalized eigenvector, y/n? "] ; 
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If [ normeigv == "y", (*then*) 

msOv = Sqrt [ Simplify [ asqr ] ]; 
sOvl = sOvl/msOv; s0v2 = s0v2/ins0v; 
asqr = 1; 

Print [ "\n(sOv(x), sOv'(x)) = " ]; 

intg = cc[sOvl] D[ sOvl, x ] + cc[s0v2] D[ s0v2, x ]; 
If [ exptrig == "y", intg = Expand[ intg, Trig -> True ] ]; 
intg = Simplify [intg] ; 
(** Print ["intg = ", intg]; **) 

If[ intg =!= 0, (*then*) (* Print [ "No, (sOv, sOv'(x)) =" ]; *) Print [ intg ]; 
intheta = InputString[ 

"Integrate (sOv(x), sOv'(x)) dx to calculate theta, y/n? "] ; 
If [ intheta == "y", (*then*) 

theta = i Integrate [ intg, x ]; 

Print ["theta = ", theta]; 

WriteString [vc , "\n theta = \n"] ; Write[vc, theta]; 
phasf = Cos [theta] + i Sin [theta]; 
sOvl = sOvl phasf; s0v2 = s0v2 phasf 
] , (*else*) Print [0] 

]; 

]; 

sOvl = Simplify [ sOvl ]; s0v2 = Simplify [ s0v2 ]; 

sOv = { sOvl, s0v2 }; 

spv = -[ - cc[s0v2], cc[sOvl] }; 

(**) 

WriteString [vc, "\n sOv = \n"] ; Write[vc, sOv] ; 
WriteString [vc, "\n spv = \n"] ; Write [vc, spv]; 
WriteString [vc , 

"\n *** sv_m = cp_m spv + c_m sOv, m = 1, 2, mmax *** \n"] ; 

(**) 

sOvlms = sOvl cc[sOvl]; s0v2ms = s0v2 cc[s0v2]; 
den = sOvlms (G22 - Qsq) + s0v2ms (Gil - Qsq) - 

cc[sOvl] s0v2 G12 - sOvl cc[s0v2] G21; 
If[ exptrig == "y", den = Expand[ den. Trig -> True ] ]; 
den = Apart [ Simplify [ den ] ]; 
WriteString [vc, "\n D = \n"] ; Write [vc, den]; 
(**) 

coef = - 2 Qsq/den; 
(**) 

Y[0] = 1; sv[0] = sOv; 
(**) 

bv[l] = i Qml D[sv[0], x] ; 

hrmthQ = InputString ["Hermiticin or Non-hermitian theory, h/n? "] ; 
WriteString [vc , 

If[ hrmthQ == "h" , " *** Hermitian ", " *** Non-hermitian "] ]; 

WriteString [vc, "theory *** \n"] ; 

simthQ = If [ hrmthQ == "h", (*then*) InputString [ 

"Simplified, Fulling or Wronskian conserving theory, s/f/w? "] , 

(*else*) "s" ] ; 
wresQ = InputString ["Write results, y/n? "] ; 

sresQ = InputString["Simplify results, y/n? "] ; 
If [ wresQ == "y" , (*then*) 
outpYm = InputString [ 

"Simple fractions. Common denominator or NO output spec, in Ym, s/c/n? "] ; 
outpsv = InputString [ 

"Simple fractions. Common denominator or NO output spec, in svm, s/c/n? "] ]; 

aprog = InputString ["Append progrzim, y/n? "] ; 

(* 
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******************* Start iterations for Y[in] and sv[in] ********************** 
*) 

For[iii =2, m <= mmxpl, m++, 

ml = m - 1; 
(**) 

cpf [ml] = coef i - s0v2, sOvl } . bv[ml]; 

If[ exptrig == "y", cpf [ml] = Expand[ cpf [ml] , Trig -> True ] ]; 
Print [ " m = " , ml ] ; 

intgnt[ml] = If [ simthQ == "s", 0, If [ OddQ[ml] && simthQ == "w", 

2 re[ cpf [ml] { cc[D[sOvl, x] ] , cc[D[s0v2, x] ] } . spv ], 
i 2 im[ cpf [ml] i cc[D[sOvl, x]], cc[D[s0v2, x]] > . spv ] ] ]; 
If[ exptrig == "y" , intgnt [ml] = Expand [ intgnt [ml] , Trig -> True ] ]; 
(**) 

If [ ml > 1, (*then*) 

For[ alpha = 1, alpha <= ml-1, alpha++, 

intgnt [ml] -= If [ simthQ == "s", 0, If [ OddQ [alpha] && simthQ == "w" , 
-1, 1 ] -[cc [sv [alpha] [[1]]] , cc[sv[alpha] [[2]]]} . D[sv[ml - alpha], 
X ] ] 

] 

]; 

If [ exptrig == "y" , intgnt [ml] = Expeind[ intgnt [ml] , Trig -> True ] ]; 
intgnt [ml] = Simplify [ intgnt [ml] ] ; 

(* 

WriteString[vc, "\n m = \n"] ; Write[vc, ml]; 

WriteString [vc , "\n integrant[m] = \n"] ; Write[vc, intgnt[ml]]; 

*) 

(**) 

cf[ml] = If[ simthQ == "s", 0, Integrate[ intgnt [ml] , x ] ]; 
sv[ml] = cpf [ml] spv + cf [ml] sOv; 

If[ exptrig == "y", sv[ml] = Expand[ sv[ml]. Trig -> True ] ]; 
(*** sv[ml] = Simplify [ sv[ml] ] ; ***) 
(**) 

Y[ml] = If [ hrmthQ == "h" , 

(*then*) ( { cc[sOvl], cc[s0v2] > . bv[ml] )/asqr, 
(*else*) (Qm2 cpf [ml] G12 asqr/(2 sOvl) + 
bv[ml] [[l]])/sOvl ] ; 
If[ exptrig == "y", Y[ml] = Expand [ Y[ml], Trig -> True ] ]; 
(* WriteString [vc, "\n Y[m] = \n"] ; Write[vc, Y[ml]]; *) 
(*** Y[ml] = Simplify [ Y[ml] ]; ***) 
(**) 

If [ m < mmxpl, (*then*) 

suml = 0; sum2 = 0; sum3 = 0; sum4 = 0; sum5 = 0; sum6 = 0; 
Do[ suml += If [ a + b + s == m, (*then*) 

Y[a] Y[b] ( sv[s] + 2 (Y[s] sv[0] - bv[s]) ), (*else*) 

], {a, 0, m - 1}, {b, 0, m - 1}, {s, 1, m - 1> ]; 
Do[ sum2 += If [ a + b + g + d + s == m, (*then*) 

Y[a] Y[b] Y[g] Y[d] sv[s], (*else*) ], 

{a, 0, m - 1}, {b, 0, m - 1}, {g, 0, m - 1}, {d, 0, m - 1}, {s, 1, m - 1} ]; 
Do[ sum3 += If[ a + b == m, (*then*) Y[a] Y[b], (*else*) ], 

{a, 1, m - 1>, {b, 1, m - 1} ] ; 
Do[ sum4 += If [ a + b + g + d == m, (*then*) 

Y[a] Y[b] Y[g] Y[d], (*else*) ], 

{a, 0, m - 1}, {b, 0, m - 1}, {g, 0, m - 1}, {d, 0, m - 1} ] ; 
Do[ sum5 +=If[a+b+g+s==m- 1, (*then*) 

Y[a] Y[b] Y[g] Qml D[sv[s], x] , (*else*) ], 

{a, 0, m - 1}, {b, 0, m - 1}, {g, 0, m - 1>, {s, 0, m - 1} ] ; 
Do[ sum6 += If [ a + b + s == m - 2, (*then*) 

Y[a] ( Y[b] ( Qm2 ( D[sv[s], {x, 2}] - Qml dQ D[sv[s] , x ] ) + 



epsO sv[s] ) - Qm2 D[Y[b], x] D[sv[s], x] - 

(1/2) qm2 ( D[Y[b], {x, 2}] - Qml dQ D[Y[b], x ] ) sv[s] ) + 
(3/4) qm2 D[Y[a], x] D[Y[b], x] sv[s], (*else*) ], 
{a, 0, m - 2}, {b, 0, m - 2}, {s, 0, m - 2} ] ; 

(**) 

bv[m] = (1/2) (suml - sum2 + ( sum3 - sum4 ) sv[0] + i 2 sum5 + sum6)] 
(**) 

]; 

t = TimeUsed [] ; 

WriteString [vc , "\n CPU time used for computation (seconds) = "] ; 

Write [vc, t - tO] ; 

(**) 

If [ wresQ == "y", (*then*) 
(* 

Simplify ajid/or writG rGsults ************************* 

*) 

For [m = 1 , m <= mmax , m++ , 

WriteString [vc , "\n m = "] ; Write[vc, m] ; 

aux = If[ sresQ == "y", (*then*) Simplify [ Y[m] ], (*else*) Y[m] ]; 
auxl = If [ outpYm == "c", Together [aux] , If [ outpYm == "s", Apart [aux] , 
aux ] ] ; 

WriteString [vc , "\n Y_m = \n"] ; Write[ vc , auxl ]; 
(**) 

aux = cpf [m] ; 

auxl = If [ exptrig == "y", Expsind[aux, Trig -> True], aux ]; 
aux2 = If [ sresQ == "y", (*then*) Simplify [auxl] , (*else*) auxl ]; 
aux3 = If [ outpsv == "c". Together [aux2] , If [ outpsv == "s". Apart [aux2] , 
aux2 ] ] ; 

WriteString [vc , "\n cp_m = \n"] ; Write[ vc , aux3 ]; 
(**) 

aux = cf [m] ; 

auxl = If [ exptrig == "y", Expand[aux, Trig -> True], aux ]; 
aux2 = If [ sresQ == "y", (*then*) Simplify [auxl] , (*else*) auxl ]; 
aux3 = If [ outpsv == "c". Together [aux2] , If [ outpsv == "s". Apart [aux2] , 
aux2 ] ] ; 

WriteString [vc , "\n c_m = \n"] ; Write[ vc , aux3 ] 
] 

]; 
(**) 

If [ aprog == "y", (*then*) << ap.mat ]; 
t = TimeUsed [] ; 
WriteString [vc , 

"\n CPU time used for computation & simplification (seconds) = "] ; 
Write [vc, t - tO] ; 

End of file Ymsvm.mat. 

The file that follows is an example of the data file for the program Ymsvm. Again it is an empty file containing c 
comments. By uncommenting appropriate pieces, one can activate input data pertaining to examples given in 
Sec. VIII. 

File Ymsvm.dat: 

Data for program Ymsvm ****) 

(**) 

(*** Example B ***) 
(* 

af = 0; 
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Rll = - ( X Cos[x]~2 + Sin[x]"2 ); 

R12 = (x - 1) Cos[x] Sin[x]; 

R21 = (x - 1) Cos[x] Sin[x]; 

R22 = - (x Sin[x]~2 + Cos[x]"2 ); 

parrepls = ■[ x -> 2 }; 

*) 

(*** Example C.l ***) 
(* 

af = 0; 

Rll = X Cos[x]~2 + Sin[x]''2; 
R12 = i (x - 1) Cos[x] Sin[x]; 
R21 = - i (x - 1) Cos[x] Sin[x]; 
R22 = X Sin[x]~2 + Cos[x]'2; 
parrepls = { x -> 2 }; 
*) 

Example C.2 ***) 

(* 

af = 0; 

Rll = - ( X Cos[x]"2 + Sin[x]"2 ); 
R12 = i (x - 1) Cos[x] Sin[x]; 
R21 = - i (x - 1) Cos[x] Sin[x]; 
R22 = - ( X Sin[x]"2 + Cos[x]"2 ); 
parrepls = { x -> 2 }; 
*) 

(*** Example C.3 ***) 
(* 

af = 0; 

Rll = - ( X Cos[x]~2 + Siii[x]'2 ); 

R12 = (1 + i)/2-(l/2) (x - 1) Cos[x] Sin[x] ; 

R21 = (1 - i)/2-(l/2) (x - 1) Cos[x] Sin[x] ; 

R22 = - (x Sin[x]~2 + Cos[x]"2 ); 

parrepls = { x -> 2 }; 

*) 

Example C.4 ***) 

(* 

af = 0; 

Rll = - ( X Cos[x]'2 + Sin[x]'2 ); 

R12 = (Cos[fi] + i Sin[fi]) (x - 1) Cos[x] Sin[x]; 

R21 = (Cos[fi] - i Sin[fi]) (x - 1) Cos[x] Sin[x]; 

R22 = - (x Sin[x]~2 + Cos[x]~2 ); 

(*fi = x;*) 

parrepls = { x -> 2 }; 

*) 

(*** Example D ***) 
(* 

af = 0; 

Rll = X Cos[x]"2 + Sin[x]"2; 

R12 = 2 i (x - 1) Cos[x] Sin[x]; 

R21 = - (1/2) i (x - 1) Cos[x] Sin[x]; 

R22 = X Sin[x]~2 + Cos[x]~2; 

parrepls = ■[ x -> 2 }; 

*) 

(*** Example E ***) 
(* 

af = 0; 

Rll = hO[x] - hl[x] ; 
R12 = h2 [x] ; 
R21 = h2[x] ; 
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R22 = hO [x] + hi [x] ; 
automatic = False; 

r[x_] := Sqrt[ hl[x]-2 + h2[x]-2 ]; 

Delta = 4 r [x] "2; 

sqrtDel = 2 r [x] ; 

signQsq = - 1; 

*) 

Excunple X 

(* 

af = 0; 

Rll = El - x~2; 
R12 = x; 
R21 = x; 

R22 = E2 - 4 x~2; 

parrepls = { x -> 3, El -> 1, E2 -> 2 }; 
*) 

End of file Ymsvm.dat. 



The file that follows contains an example of a program that can be appended to Ymsvm by answering y to the 
dialog question Append program, y/n?. This program is pertinent to Example E in the data file Ymsvm.dat. It 
computes and prints both the formulas in Eq. (174) and all numerical results given in TABLE I in [l|. For each of two 
eigenvalues (for the minus or plus sign in Eq. (123)), the program Ymsvm should be run twice for two possible answers 
to the dialog question Normalized eigenvector, y/n?. In both cases one should take mmax = 2 and answer n 
to the dialog questions Write results, y/n? and Simplify results, y/n?. Otherwise, the program would 
spend very long time in an attempt to simplify Y2{x) and c^(a;) which in any case are too complicated to be presented. 



File ap . mat : 

(***** Appended computation for program Ymsvm, Example E *****) 
(**) 

aux = Together [ Y[l] ]; 

WriteString[ vc , "\n Y_l = \n"] ; Write[ vc, aux ]; 
(**) 

aux = Together [ cpf[l] ]; 

WriteString [vc , "\n cp_l = \n"] ; Write[ vc , aux ]; 
(**) 

dO[x_] := 1/(4 x-2) + 4/x-4 + 38/x~6 + 748/x-8; 
dl[x_] := l/x-2 + 2/X-4 + 19/x-6 + 374/x-8; 
hO[x_] := -1 - k'2 + dO [x] ; hi [x_] := 2 (om + l/x~2) 
(**) 

parrepls = { x -> 55, k -> 4/100, om -> 26041/10~7 >; 
WriteString[ vc , "\n"] ; Write[ vc , parrepls ]; 

"\n Q = "] ; Write[ vc , N[ Q /. parrepls ] ]; 
"\n epsO/2 = "] ; Write[ vc , N[ epsO/2 /. parrepls ] ]; 
"\n Y_l = "] ; Write[ vc, N[ Y[l] /. parrepls ] ]; 
"\n Y_2 = "] ; Write[ vc, N[ Y[2] /. parrepls ] ]; 
"\n cp_l = "] ; Write[ vc, N[ cpf[l] /. parrepls ] ]; 
"\n cp_2 = "] ; Write[ vc, N[ cpf [2] /. parrepls ] ]; 



h2[x_] := -1 + dl[x] ; 



WriteString [ vc , 
WriteString [ vc , 
WriteString [ vc , 
WriteString [ vc , 
WriteString [ vc , 
WriteString [ vc , 



End of file ap . mat . 



[1] A. A. Skorupski, Phase Integral Approximation for coupled ODEs of the Schrodinger type, larXiv:0710.5868 [math-ph]. 



